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A particular mode of isotachophoresis (ITP) employs a pressure-driven flow opposite to the sam- 
ple electromigration direction in order to anchor a sample zone at a specific position along a 
channel or capillary. We investigate this situation using a two-dimensional finite-volume model 
based on the Nernst-Planck equation. The imposed Poiseuille flow profile leads to a signifi- 
cant dispersion of the sample zone. This effect is detrimental for the resolution in analytical 
applications of ITP. We investigate the impact of convective dispersion, characterized by the 
area-averaged width of a sample zone, for various values of the sample Peclet-number, as well as 
the relative mobilities of the sample and the adjacent electrolytes. A one-dimensional model for 
the area-averaged concentrations based on a Taylor- Aris-type effective axial diffusivity is shown 
to yield good agreement with the finite-volume calculations. This justifies the use of such simple 
models and opens the door for the rapid simulation of ITP protocols with Poiseuille counterflow. 

I. INTRODUCTION 

Isotachophoresis (ITP) is a special mode of electrophoretic transport that has already 
been described by KohlrauscbP. It has found widespread applications in the analytical 
sciences (for an overview, see reference^) and is mainly used as a technique for the pre- 
concentration and separation of samples. ITP relies on consecutively stacking a high 
mobility leading electrolyte (LE) and a low mobility trailing electrolyte (TE) in a cap- 
illary or channel as sketched in figure [T] Upon application of an electric field the ions 
arrange in the order of their electrophoretic mobility, forming a sharp transition zone 
between them that migrates along the capillary according to the electrophoretic velocity 
of the ions. Accordingly, there is a corresponding change in electric field across this tran- 
sition zone, such that the high and low mobility species move at the same speed. A high 
mobility ion diffusing backwards across the transition zone into the low mobility region 
will experience this higher electric field and is thus transported back to its own region 
and vice versa for a low mobility ion diffusing into the region of high mobility ions. It is 
this balance between diffusion and electrophoretic migration that determines the width 
of the transition zone. Depending on whether the stacking occurs for the anions or the 
cations, anionic or cationic ITP is observed. Sample ions with an electrophoretic mobility 
in between the LE and TE ion mobilities are sandwiched between the two electrolytes 
and can form a zone of only a few micrometers width, depending on the total amount of 
sample. 

Very early in the development of electrophoretic separation techniques researchers ap- 
plied a counterflow opposing the electrophoretic sample migration to increase the time 
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span available for the separation of analytes or to reduce the length of the separation 
passage^ . While this principle was first applied to capillary electrophoresis, Everaerts et 
alP introduced a counterflow to balance the sample migration in ITP experiments. Orig- 
inally, pressure-driven counterflow was applied, but later electroosmotic counterflow was 
considered too^. In the past decades, counterflow-balanced ITP has found widespread 
applications, for example in the context of sample preconcentration for improving the 
sensitivity of capillary electrophoresis^ - ^ or to control the elution of analytes^. 

Despite the importance of the method, it is rather poorly understood how the coun- 
terflow affects the analyte distribution in a sample focused by ITP. Urbanek et al.^ 
have noted in their experiments that pressure-driven counterflow increases sample dis- 
persion. A few efforts have been made to compute the effects of counterflow on ITP. 
The one-dimensional models that have been develope d 11 ! * 12 ! may be suitable for flow pro- 
files resembling a plug flow (such as electroosmotic flow in capillaries that are much wider 
than the Debye layer thickness), but do certainly not account for the complex convection- 
diffusion-electromigration phenomena occurring when a Poiseuille counterflow is applied. 
More successful one- dimensional models considering convective dispersion in ITP were 
developed with the aim to describe dispersion in a setting with a mismatch in electroos- 
motic flow™ 21 . 

In this article, we present the first computational fluid dynamics (CFD) study of sample 
dispersion occurring when ITP is balanced by a Poiseuille counterflow. For this purpose 
we have numerically solved the coupled Nernst-Planck and charge conservation equations 
using a finite-volume method. We analyze the most important dependencies of the width 
of the area-averaged sample distribution on the dimensionless parameters governing the 
problem. We show that in ITP with short sample zones and at large enough (but not 
unrealistic) Peclet numbers the flow considerably broadens the sample distribution com- 
pared to pure ITP without counterflow. Furthermore, we show that a one-dimensional 
Taylor dispersion model reproduces the data obtained with CFD reasonably well. 



II. GOVERNING EQUATIONS 



The migration of charged species in an electrolyte under an external electric field is 
governed by convection, diffusion and electromigration. The mass conservation of the 
ionic species leads to the Nernst-Planck equation, 

dC 

+ v • [(v + z^e) d - Avcy = o. (i) 

Here Cj, Di, fi i: Zi denote, respectively, the concentration, diffusion coefficient, mobility 
and valency of the i th ionic species. Our model system comprises four ionic species: ions in 
the TE, sample and LE having a valency of the same sign and a common counter-ion with 
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valency of opposite sign. The subscript i denoting these species is chosen as t, s, I, and 
0, respectively. The mobility and diffusivity are related via the Nernst-Einstein relation 
fii = DiF/RT, where F is the Faraday constant, R the gas constant and T the absolute 
temperature. V is the flow velocity field and E = — V0 the electric field. The electric 
current density and charge density are defined as j = FJ2i z i^i an d p e = FJ2i z iCi, 
where Nj = — DjVCj + Cj(V + z^E) denotes net flux of the i th ionic species with Zi — 1 
for i = t, s or I and zq = —1. 
The charge conservation equation is 

t + V.i = 0. (2) 

Under electro-neutrality^, i.e. p e = 0, the time-dependent and convective terms in (|2| 
vanish. By multiplying equation by Fzt and taking the sum over all species the 
reduced charge conservation equation can then be written as 

V>E)=V-(F^A^VC<), (3) 



where the ionic conductivity is v — (F^^izfCi). The right-hand side term in ([3]) is 
known as the diffusion current and its contribution is insignificant at all locations, except 
for the transition zones where large concentration gradients occur. For microchannels 
with charged walls, electroneutrality is violated within a region of the order of the Debye 
screening length, Xo, much smaller than the channel dimensions. In such cases equation 
(|3| is retained in the channel bulk, while the effect of the charge layer is incorporated 
via an effective wall-velocity boundary condition taking into account the electroosmotic 
flowff^. Since here we focus on sample dispersion by Poiseuille flow we assume uncharged 
walls. 

In this article we consider a two-dimensional situation, i.e. the isotachophoretic trans- 
port and the Poiseuille flow occur between parallel plates. Translated to realistic mi- 
crofluidic setups this means that very shallow channels with a width much larger than 
their depth are assumed. We impose a Poiseuille counterflow with an average speed equal 
and opposite to ITP migration, as 

u(Y) = -6U< TP l (l-l), (4) 



H \ H / 

where H is the depth of the channel and U ITP is the velocity of ITP transport for the case 
that no counterflow is applied. The situation is depicted in figure [TJ Thus, on average 
the electromigration of sample ions is just compensated by convection, meaning that the 
average speed of the ITP transition zone (the "interface" between LE and TE) will be 
zero. 

We non-dimensionalize the X- and F-coordinates by the depth H of the channel, i.e. 
set x = X/H and y = Y/H. Further, all concentrations are non-dimensionalized by the 
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FIG. 1. Sketch of the situation considered. A single sample zone is sandwiched between a leading (LE) and 
trailing (TE) electrolyte. The sample ions migrate from left to right at velocity U ITP in an applied electric field. 
A pressure driven flow from right to left is applied that exactly counters the electromigration, such that the 
sample plug is stationary. The channel length, L = L t + L s + Li, is the sum of lengths occupied by the leading 
and trailing electrolytes as well as the sample. We will consider L s C Li, Li. 



bulk LE concentration, Cf°, and the electric potential by <fto = RT/F. We choose the 
isotachophoretic velocity U ITP as velocity scale. Time is scaled by r = H/U ITP , i.e. 
i = t/r. The Nernst-Planck equation ([!]) can then be written in non-dimensionalized 
form as 



dci dd Di 1 Di 1 2 



0. 



(5) 



where & 



Cj/Cj 00 are the non-dimensional concentrations and the Peclet number is 



defined as Pe = U ITP H/D S . The charge conservation equation ^ becomes 



Q + ^ , ^ C s + — — ; — —C t V(p 



A + A) 

ziDi + zqD 
Di + Do 



A + D 
V 2 r 



z s D s + zqD q z t D t + z D 
°i H 7=r c s H ^— ^-Q 



(6) 



z x Di + z D ziDi + z D 
Under no counterflow, all the zones move with a constant isotachophoretic velocity 
jjitp _ ^ t E t — jji s E s = fiiEi, where Ei (i — t,s or /) are the local electric field strengths 
in the TE, sample or LE zone, respectively, given as 



Ei = En 



(7) 



where U = L^/L. L,i (i — t,s or I) is the portion of the channel filled by the respective 
electrolyte and L is the total length of the channel. E = $/L is the average electric 
field due to the voltage drop $ along the channel. Here we have assumed that the sample 
forms a distinct zone where its concentration remains constant and the overlap region 
between the zones are negligible. Since in ITP applications the channel is long compared 
to the section occupied by sample we take the limit l s — > in equation Q, which then is 
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also valid for the electric field in the TE and LE even in the case of a dispersed sample 
zone. Also, we assume that the sample zone is situated in the middle of the channel, i.e. 
we choose k = l t = 0.5. 

From the conservation of electric current along with electro-neutrality and the negligible 
diffusive flux at the far-field, a relationship between the TE and LE concentrations can 
be obtained as 

c t °° = m + fig fit gf = m + fio fi a 

C7° A** + ^0 M C?° + /^o M 
where the suffix 00 stands for the concentrations far away from the interfaces between 
the electrolyte j 1 ^ . We will refer to these relations as the Kohlrausch conditions. Note 
that the second of these relations is only valid in the case where the sample concentration 
forms a distinct zone with constant concentration, i.e. when it is only marginally affected 
by diffusive and convective dispersion. 

In solving the Nernst-Planck equation, the net ion fluxes through the channel walls are 
set to zero, i.e. Nj • n = for i = t, I or s, where n is unit outward normal. The electric 
potential is subjected to insulating boundary condition (V0- n = 0) along the wall. Both 
the left and right boundaries of the computational domain are placed sufficiently far 
away from the sample zone such that the distribution of sample ions, c s , is not affected 
by their presence. The concentrations far away from the transition zones are governed 
by the Kohlrausch condition ([8]). 

Due to the presence of different electrolytes in different portions of the channel, a sharp 
change in conductivity will appear across the transition zones. If there is no bulk flow, 
the width of such a transition zone is of the order oP^ 

A* = % ,9) 

where AE is the change in electric field strength across the zone. Upon application of a 
voltage drop $ across the channel, a locally uniform electric field appears sufficiently far 
away from each transition zone. This allows for calculation of the electric potential across 
the channel via charge conservation ([3]). In other words, owing to the electro-neutrality 
assumption we are spared the effort of solving the Poisson equation from which the electric 
potential is usually derived. 

The electroneutrality assumption^ relies on the fact that its inevitable violation dic- 
tated by a changing electric field due to Gauss's law, p e = e e divE, is much smaller than 
the total amount of charges present, ~ F|zj|cj. For ITP applications the ratio between 
the two is typically of the order of 10 -4 — 10~ 5 (cf. supplementary material fo j 16 * 21 *) which 
allows us to safely neglect the impact of p e j F on the concentrations. However, this is not 
the only place where the charge density enters, since in terms of Maxwell stresses it also 
contributes as a body force, ~ p e E, in the momentum equation for the flui d 16 * 22 *. This 
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issue of Maxwell stresses in the context of ITP in a similar setting was discussed irpH, 
where it is concluded that the body force term becomes important for large applied fields 
with very narrow transition zones and hence large field gradients, something that is not 
achieved with the strong convective dispersion considered here. 

In order to make quantitative predictions using this model, the electrophoretic mobili- 
ties of all species, concentration of the leading electrolyte, voltage drop across the channel, 
and depth of the channel must be specified. Note that the TE concentration can be eval- 
uated by using the Kohlrausch conditions (j8|). The electric field strengths at the ends of 
the channel can be adjusted via relation (|7|). We use e e = 78.5 x 8.85 • 10~ 12 C/(Vm), 
T = 300 K and O = 25.9 mV. 

III. NUMERICAL METHODS 

We developed a computer code to compute the governing equations based on the numer- 
ical algorithm as described below. The non-dimensional equations for ion transport and 
electric field are computed in a coupled manner through the finite volume method^. The 
computational domain is subdivided into a number of control volumes. When the elec- 
tromigration in the Nernst-Planck equations dominates the electro-diffusional transport 
of ions, the transport equations show hyperbolic characteristics. Due to the hyperbolic 
nature of the ion transport equations, we adopt a higher-order upwind scheme QUICK^ 
(Quadratic Upwind Interpolation Convection Kinematics) to discretize the electromi- 
gration and convection terms in the ion transport equations. The diffusion flux at the 
control volume interfaces is estimated by a linear interpolation of variables between the 
two neighbors to either sides of the control volume interfaces. Details of the spatial dis- 
cretisation scheme can be found in the supplementary material. An implicit first-order 
scheme is used to discretize the time derivatives present. At every time level we solve 
the ion transport equations iteratively, as the ion transport equations and the charge 
conservation equation are coupled. The iteration procedure starts with a guess for the 
electric potential at each cell center. At every iteration, the elliptic PDE for the charge 
conservation equation is integrated over each control volume through the finite volume 
method. The discretized equations are solved by a line-by-line iterative method along 
with the successive-over-relaxation (SOR) technique. The iterations are continued until 
the absolute difference between two successive iterations becomes smaller than the tol- 
erance limit 10~ 6 for concentration as well as for the electric potential. The details on 
discretization of the governing equations are provided in the supplementary material. 

A steady state solution is achieved by taking a sufficient number of time steps until the 
concentration distributions remain unchanged with time. The initial condition for ITP 
with counterflow is governed by the solution of the corresponding ideal ITP case (without 
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convection). The solution for ideal ITP is also obtained based on the algorithm as de- 
scribed above and the parabolic velocity profile is imposed after a steady isotachophoretic 
zone is formed. We find that the steady-state state is archived after a short transition 



phase. The time evolution and formation of steady-state is illustrated in figure S4 of the 
supplementary material. However, in this paper all results presented correspond to the 
situation where a steady state has been reached. 

As the variables within the transition zone between two electrolytes change more rapidly 
than elsewhere, a nonuniform grid spacing along the x-direction and uniform spacing 
along y-direction is chosen. The grid size Sx within the transition zones is relatively 
small and is increased with a constant increment as we move away from the transition 
zone. The smallest Sx is chosen such that the Courant-Friedrich-Levy (CFL) criterion is 
satisfied. For the present computations, Sx is taken to be 0.005 for the finest grid and it 
is 0.01 for coarse grids. Here Sy is taken to be 0.005 and St = 0.001. A grid independency 
test and validation of our algorithm by comparing with the analytical solution for ideal 



ITP (Goet et al.^j) is presented in figure S3 in the supplementary material. The efficiency 
of the present numerical code in resolving the transient zones for peak and plateau-mode 



ideal ITP in steady-state is illustrated in figure S5 of the supplementary material. 



IV. TAYLOR-ARIS DISPERSION MODEL 



To analyze the dispersion effect on the sample when a counterflow is applied, we consider 
a one-dimensional model based on Taylor-Aris dispersiorpSES Taylor's original work was 
concerned with dispersion solely due to the pressure driven flow through a channel. In 
that case the mechanism at work is the interplay of convective dispersion and lateral 
diffusion. Convective dispersion alone would result in a dispersion linear with time, but 
its effect is limited by the lateral diffusion over the cross section of the channel. Combined, 
these phenomena result in a dispersion proportional to the square root of time, such that 
the area-averaged sample distribution can be considered to spread under the influence of 
an effective axial diffusion. For two reasons it is a priori unclear whether such a Taylor- 
Aris dispersion model can be applied to the situation considered in this article. First, 
electromigration appears as an additional transport process that results in a sharpening 
of the sample distribution. Secondly, an assumption implicit to the Taylor-Aris model 
is the infinite time limit, i.e. the model is valid only on time scales large compared 
to the time a molecule takes to diffusively sample the channel cross section. In turn, 
this means that it is usually applicable only to comparatively broad (in axial direction) 
sample distributions with a/H 3> Pe, where a is the length scale of a transition zone. For 
this reason one would a priori not expect that a description based on an effective axial 
diffusivity is applicable. Nevertheless, such models have been successfully used previously 
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in similar situation d 13 H 15 N 16 [ Therefore we will investigate the applicability of this model 
to the case considered here. 

The core of the Taylor-Aris dispersion model is the effective axial diffusion due to 
convection. In particular, the diffusion coefficient is replaced by a term dependent on the 
Peclet number squared 

where i = t, s or I. For a parallel-plates channel the value of the constant /3 is 1/210. 
The dispersion model is based on the area-averaged ion concentrations 
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Likewise, we write area-averaged transport equations for the ion concentrations in which 
the diffusivity is replaced by the effective axial diffusivity according to Taylor and Aris 
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The axial electric field is determined by solving the corresponding one-dimensional charge 
conservation equation assuming electro-neutrality throughout 

d 



dX 



(ve x ) = fY j dJ^ 



(13) 



where v is the area-averaged conductivity. The equation for the electric field is coupled 



with the transport equations (12). The numerical solution of the species transport equa- 
tions above is achieved by employing a higher-order accurate upwind differencing scheme. 
The equation for the electric field is solved along with the ion transport equations in a 
coupled manner as outlined before. 



V. RESULTS AND DISCUSSION 



We consider the LE and common ion diffusivities as fixed values throughout this study, 
setting Di = Dq = 7.0 x 10~ 10 m 2 /s. The diffusivities of the TE and sample ions are 
prescribed by specifying the two non-dimensional parameters k\ = Di/D t and k 2 = 
Di/D s . The corresponding mobilities are related to the diffusivities via the Nernst- 
Einstein relation, ^ = DiF/RT. Since the sample mobility lies between the mobilities 
of TE and LE, these parameters can be varied subject to the restriction k\ > k 2 > 1. 
The scale of the applied electric field is set to Eq = 10 5 V/m from which the field in 
the respective buffers can be inferred using equation ([7]). The bulk LE concentration is 
taken as = 10~ 3 M , so that the Debye layer thickness (~ 10 nm) is much smaller than 
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the depth of the channel considered here (> 5/xm). The bulk TE concentration can be 
inferred from the Kohlrausch condition rf8J), which in our example becomes C^/Cf = 
2/(1 + ki). A corresponding upper bound can be obtained for the maximum of the 
sample concentration, which due to dispersion will usually be much smaller in the cases 
considered here, however. 

In ITP experiments one considers two modes of operation, peak and plateau mode, 
depending on the amount of sample present in the system. "Plateau mode" refers to the 
case where the sample-zone is wide compared to the transition zones, i.e. the sample 
concentration distribution forms a constant plateau with blurred boundaries towards LE 
and TE. "Peak mode" refers to the other extreme of a very short sample zone, where 
the two transition zones at both sides of the sample overlap or more precisely, when the 
sample is entirely within the diffused interface between LE and TE. In this situation the 
influence of the analytes on the overall conductivity is typically negligible. Everything 
else being equal, the difference between peak and plateau mode lies in the amount of 
sample present. In the present case, owing to the strong additional dispersion due to 
convection, a larger amount of sample will still lead to a pronounced peak instead of a 
plateau, as will be elaborated below. In order to distinguish this regime from peak mode 
as observed without convective dispersion, we will call this situation "dispersed plateau 
mode". This is the regime that we will mainly, but not exclusively, consider presently. 
It is different from plateau mode in that the area averaged sample distribution shows 
a pronounced peak and differs from peak mode in that the sample strongly affects the 
overall conductivity. 

In order to obtain comparable results for the distortion of the sample region due to the 
applied Poiseuille flow, we fix the area-averaged amount of sample 

/oo 
dXC s (X) (14) 
-oo 

present in the channel. In that way the area-averaged distribution without convective dis- 
persion is independent of the channel height. Note that the exact form of the distribution 
will depend on the choice of other parameters such as k\ and ki due to the different con- 
centrations dictated by the Kohlrausch conditions QSJ) or the width of transition zones Q. 
Unless stated otherwise we will mainly use the values of C s> o = 40 /zmol/m 2 = 40 /zm • Cf° 
and C S fl = 4/im • Cf° throughout; the former value leads to a distinct plateau without 
convective dispersion while it represents dispersed plateau mode for wide enough chan- 
nels, while the latter corresponds to a sample distribution in the form of a Gaussian peak 



without convective dispersion, cf. figure [S5] in the supplementary material. 

We remark that without convective dispersion the width of the sample zone in plateau 
mode can easily be inferred from the Kohlrausch condition ^ and the amount of sample 
present, neglecting the width ^ of the diffusive zones, to L s ps C s /C^° = C s /Cf° ■ (1 + 
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FIG. 2. Sample concentration obtained in the 2D model for fci (= Hi//J,t = Di/Dt) = 3, k2{= Hi/lAs = Di/D„) = 2 
with channel depths of (a) H = 10 /im (Pe = 40) and (b) H = 30 fj,m (Pe = 120). The amount of sample is 
Cs/Cj" — 40 /im. Both coordinates axes show the length in units of 10~ 5 m, but were rescaled differently for 
better visualization of the sample distribution. 
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FIG. 3. Effect of the channel height on the area-averaged sample concentration profile and comparison with the 
corresponding profile obtained in the Taylor- Aris dispersion model for fixed values of the diffusivity ratios i.e., 
ki(= fii/ fit — Di/Dt) = 3 and &2(= fil/fis = D i/D s ) = 2 when the depth of the channel is chosen as H — 10 /im 
and H = 30 ^m, respectively. The amount of sample is C s /Cf° = 40 fim. 



k-ij/2, which gives a baseline for the width of the sample distribution. 

In figure [2] the corresponding sample concentration profiles are shown for two different 
channel heights. As mentioned, all results correspond to a situation where a steady state 
has been reached. Figure [3] shows the corresponding area-averaged sample concentra- 
tions, both derived from the 2D model and from the ID Taylor-Aris dispersion model 



(12) with effective diffusivity (see supporting information for a combined graph of the 
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TE, sample and LE distributions in this situation). From these figures it becomes appar- 
ent that whether a plateau or peak is obtained now becomes dependent on the channel 
height. For the shallow channel a plateau is observable in the area-averaged concentra- 
tions, contrasting the distinct peak observed for the wider channel, characteristic for the 
"dispersed plateau mode". Since in both cases the ITP velocity is the same, the larger 
dispersion is a direct consequence of the larger Peclet number for the wider channel. We 
also note that the ID model based on the effective Taylor-Aris diffusion captures the 
area-averaged distribution quite well. This is particularly true for the narrower channel. 
In the wider channel the larger discrepancy is due to the slight focusing of the sample 
in the center region of the channel, which effectively results in a narrower peak than 
determined from the ID model. We will further investigate this below by varying the 
total sample concentration present in the channel. 

Before doing so, we investigate the area-averaged distributions obtained with the two 
models for different ratios of the diffusivities, cf. figure [4] (a) when the sample amount 
is C s /C;°° = 40 /xm. Corresponding results for the reduced sample amount i.e., C s /Cf° = 
4 /im is shown in figure [4] (b). In the cases where the sample diffusivity is very close to the 
diffusivity for either LE or TE one expects a very wide transition zone between the sample 
and the respective electrolyte, since hardly any electromigrative sharpening is present. 
On the other hand, for the transition towards the electrolyte with markedly different 
diffusivity a much narrower zone is expected. This is indeed seen in figure [1] (a) and (b). 
We remark that the three situations depicted correspond to different Peclet numbers, 
since the parameter varied is the sample diffusivity. Nevertheless, the total width of 
the sample zone is affected more strongly by the long tails towards the electrolyte with 
similar diffusivity than by the Peclet number. We will further investigate this fact below. 
Also note that again the Taylor-Aris model quite successfully captures the area-averaged 
concentration profiles, in particular towards the electrolyte with similar diffusivity. The 
sharper profile towards the electrolyte with dissimilar diffusivity in the 2D compared with 
the ID case is again due to the more detailed 2D structure with a slight focusing in the 
channel center emerging in this model. Comparing figure [4] (a) and (b) we find that in the 
dispersed plateau mode the profiles for the sample distribution are largely independent 
of its amount. 

To further investigate the validity of the ID model, the total amount of sample present 
in the channel is varied. Figure [5] shows the 2D profiles obtained for amounts ranging from 
1/2 to 5 times C Si o, the amount considered so far. For the smallest amount of sample, one 
again observes the focusing of the sample in the central region of the channel. Contrasting 
that, for the largest amount of sample present, the situation may be considered as plateau 
mode. The corresponding area-averaged concentration profiles are shown in figure [6j both 
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FIG. 4. Area-averaged sample concentration profile obtained in the 2D simulation for different values of the 
diffusivity of sample electrolyte i.e., fe(= A*;/Ms — Di/D 3 ) — 1.1 (Pe = 55), &2 = 2 (Pe = 100), = 2.7 
(Pe = 135). Here the LE to TE diffusivity ratio i.e, ki(— ^i/fJ.t = Di/D t ) = 3 and channel depth is H — 25 fim. 
The amount of sample is (a) C s /C;°° = 40 /jm; (b) C s /C;°° = 4^m. A comparison with the corresponding profile 
obtained in the Taylor- Aris dispersion model is also made. Dashed lines represents the results obtained from the 
2D model and the solid line represents corresponding results from the ID model. 

for the 2D and ID cases. Again, the ID model is surprisingly effective in predicting the 
profiles, with the quality of the predictions deteriorating for small sample amounts due 
to the strong 2D structure of the profile, leading to a sharper focusing than predicted by 
the ID model. 

In order to quantify the dispersion, we define the width of the sample zone as the second 
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FIG. 5. Sample concentration obtained in the 2D model for fci(= fti/l-it = Di/Dt) = 3 and fc2(= i-ii/fis = 
Di/Ds) = 2 when (a) C s = 5C s ,o, (b) C s = 2C S>0 , (c) C s = lC s ,o and (d) C s = 0.5C S , with C sfi /Cf° = 40 /im. 
The channel depth is H = 25 /im. Note that both the coordinate axes are multiplied by a factor 10 _5 m. 
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where the center of the distribution is defined as 
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FIG. 6. Comparison of the area-averaged sample concentration profiles obtained in the 2D simulation with the 
corresponding profile obtained in the Taylor- Aris dispersion model for fixed values of /i;//i t = Di/D t ) = 2 
and fc 2 (= Mi/Ms = D t /D a ) = 1.5 when (a) C s = 5C S , , (b) C s = 2C Sj0 , (c) C 3 = 1C S , and (d) C s = 0.5C S , with 
C s ,o/C^° = 40 /im. The depth of the channel is H — 25 fim. The solid line represents the results obtained from 
ID model and dashed lines from 2D model. 

Additionally, we also measure the skewness of the distribution through the third nor- 
malised moment as 



to further quantify the similarity between results for the sample distribution obtained 
within the 2D model and the ID Taylor-Aris dispersion model. The value of skewness 
may be positive, negative or zero depending on whether the sample electrolyte is skewed 
to the right (positive skew), to the left (negative skew) or symmetric (zero skew). 

The stage is now set for a more systematic study of the influence of the individual 
parameters on the form of the sample distribution. We use the width of the sample zone, 
ex, and its skewness, 71, in figures |HJ [7] and [9] to investigate the dependence of the dispersion 
on diffusivity ratios of the individual electrolytes, the amount of sample present and the 
channel depth, respectively. 

In particular, figure [7] shows the distribution's width and skewness for several values of 
ki when k 2 is varied in the range 1 < k% < k±. In order to compare the results we have 
plotted them using the rescaled variable (k 2 — l)/(ki — 1) so that the curves span the 
same range in parameter space. In this figure the remark already made in the discussion 
of figure [4] is quantified: the sample dispersion is largest in the cases where the sample 
mobility is close to the mobility of one of the other electrolytes. This is reflected both in 
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FIG. 7. Variation of the (a) standard deviation and (b) skewness of the sample distribution with &2(= Mi/Ms = 
Di/D 3 ) for a fixed value of M'/m* — Di/D t ) = 2, 3, 4 and 5. The amount of sample is Cs/Cf° = 40 /im. 

Lines with unfilled symbols represent the results from the ID model (Taylor- Aris model), lines with filled symbols 
those from the 2D model. 



the increase of width, <r, as well as its skewness, 71, for values of k 2 close to either of the 
extreme cases 1 or k\. The sign of the skewness reflects the fact that for k 2 close to 1, 
i.e. the diffusivity is very close to the LE diffusivity, the sample distribution has a long 
tail reaching into the LE, and similarly when the diffusivity is close to the TE diffusivity, 
i.e. k 2 close to fci, the tail reaches to the left into the LE. We find the skewness to 
be almost zero, i.e. a symmetric sample distribution, when its mobility, fi s , is close to 
the harmonic average of [it and We note that the individual curves corresponding 
to different values of fci almost coincide showing that the dependence on k\ is relatively 
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FIG. 8. Variation of the (a) standard deviation and (b) skewness of the sample distribution with fca(= fj,i/fj, s = 
Di/D 3 ) for a fixed value of ki(= m/fJ,t = Di/Dt) = 3. The Peclet number is Pe = 50&2, the depth of the channel 
H — 25 fim. The amount of sample is C s /C;°° = 40 fim and 4^m. Lines with unfilled symbols represent the 
results from the ID model (Taylor- Aris model), lines with filled symbols those from the 2D model. 



weak when considering the rescaled k-i- values. 

To further investigate the universality of the distribution width and skewness curves, 
we show in figure [8] the corresponding curves for different amounts of sample present in 
the system. In particular, we show results for a sample amount of C s /Cf° = 40 /im and 
C s /Ci° = 4 /im and find that both the width as well as the skewness of the distributions 
do not significantly change. As noted earlier, these two situations corresponds to plateau 
and peak mode in a situation of ITP without additional convective dispersion. This 
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FIG. 9. Variation of the standard deviation of the sample distribution with the channel depth, where k 1 (= 
H>l/lM = Di/Dt) — 3, fe(= Mi/Ms — Di/Ds) = 2 and Cs/Cf = 40/im. The Peclet number varies with the channel 
depth as Pe = 100 • H/(25/Ltm). 



points towards the fact that the form of the distribution in dispersed plateau mode is not 
strongly dependent on the amount of sample present, the difference being mainly a scale 
factor proportional to the total amount, as already hinted at in the distributions shown 
in figure |4| 

Both in figure [7] and [8] we show data obtained within the 2D calculation as well as for 
the effective ID Taylor-Aris dispersion model. As can be seen, there is good agreement 
between the two calculations corroborating the use of the simplified model for predicting 
area-averaged sample distributions in cases with convective dispersion. Nevertheless, it 
is important to keep in mind that the actual sample distribution within the channel will 
usually be richer in features in the 2D case, in particular we mention the focusing of the 
sample in the channel center observed for example in figures [2] and |5j 

Another interesting aspect is how the width of the sample zone changes when varying 
the channel depth H, and hence the Peclet number, as displayed in figure [9j From the 



effective diffusivity, eq. (10), one expects a sample zone width that scales as ~ Pe 2 in 
dispersed plateau mode. For large enough H, such that the dispersed plateau mode is 
indeed observed, this is reproduced well by the ID model, as it should be. However, in 
the region shown this limit is not fully reached yet and only for values of H larger than 
30 /im the points follow a ~ H 2 for the ID results. As can be seen from the figure, 
the 2D model predicts less sample dispersion, a fact that can again be attributed to the 
focusing of sample in the channel center. Figure M (b) shows that the skewness of the 



distribution is correspondingly amplified by the increasing channel height. 
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VI. CONCLUSIONS 

We have investigated sample dispersion in situations in which electromigration in iso- 
tachophoresis is balanced by a Poiseuille counterflow. The focus was put on small sample 
amounts, i.e. situations when the area averaged sample distribution shows a peak. It is 
found that at large Peclet numbers, the sample dispersion is significantly increased com- 
pared to a situation without counterflow. Since the amount of sample is such that without 
the counterflow present one would typically obtain a short sample zone in plateau mode, 
we have termed this regime "dispersed plateau mode" . One of the findings of this study is 
summarized in figures [7] and [8] which show the dispersion of the sample zone as a function 
of the sample ion mobility. In particular, the diagram shows in which window the mobility 
of the sample can be chosen relative to the mobility of the surrounding electrolytes before 
significant broadening of the sample zone is observed. In particular, we find that the area 
averaged sample distributions are only weakly dependent on k\ = [iij \it but vary strongly 
with k 2 = Willis ( or more precisely with the rescaled value (k 2 — l)/(&i — 1))- These re- 
sults are also observed to be relatively insensitive to the amount of sample present in the 
system, cf. figure [8j with the general form of the distribution remaining the same but 
the peak-height scaling with the total amount. These curves can thus be used to asses 
the degree of dispersion and skewness expected in an experiment. Similar results hold 
for the case in which it is desired to focus two sample zones stacked behind each other 
without the risk of too strong intermixing. 

Another key result of this study is that a ID area-averaged model for the sample dis- 
tribution based on Taylor- Aris dispersion agrees well with the more detailed 2D model 
in many situations. In view of the original derivation of the effective diffusion coefficient, 
the quality of the agreement is surprising, since the complex interplay of convection, 
electromigration and diffusion that shapes a transition zone in ITP is not accounted for. 
Since the computational cost of the ID model is much less than that of the 2D model, this 
opens the door for fast simulations of ITP processes that are superimposed by convection. 
This could be of considerable relevance for a number of processes involving sample pre- 
concentration and separation. However, it is important to stress that the detailed spatial 
distribution of the sample ions generally is more complex than the simplified ID picture 
suggests. In particular we observe in the 2D calculations that the sample has a tendency 
to become concentrated in the channel center. This will have an impact in situations 
where different ionic species are made to react within an ITP experiment, as for example 
described in Goet et al.^ and Bercovici et alP^, when long reaction times are attained 
via a Poiseulle counterflow. 
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SUPPLEMENTARY MATERIAL 

SI DISCRETISATION OF GOVERNING EQUATIONS 

The Nernst-Planck equation can be written as 



where F 



^ D s Pe dx ' ^ 



-z. 



Di 1 



-V 2 c = (18) 
and c is the concentration of the i th 



' D s Pe dy 

ionic species. The computational domain is sub-divided into a number of elementary 
rectangular cells Qp with area dVtp whose sides are dx p and dy p . The ion transport 



equations, when integrated over a cell ilp (cf. figure SI), yields the discretised form to 
advance the solution from k th time step to (k + l) th time step 
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(19) 



Here n, s, e and w refer to the northern, southern, eastern, western face of the cell (cf. 



figure SI ). An implicit first-order scheme is used to discretise the time derivatives present. 
The electromigration and convection terms are discretised through a higher order upwind, 
QUICK scheme as follows 



(F e c e - F w c w )dy P + (G n c n - G s c s )dx P 
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dy P 
dxp 
dxp 

(20) 

where the operator [[a, b}} yields the larger of a and b. The diffusion flux at interfaces 'e' 
and 'u;' are evaluated as 
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-c N + -^cp - -c 3 )[[G n ,0]] - (-c N + -c P - -c NN )[[-G n ,0)) 



and 
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(21) 



(22) 



0.5(dxp + dxw) 

A similar procedure is adopted for estimating the variables at the other cell faces £ n' 
and V. Note that the big letter subscripts denote the cell centers in which variables 
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F > F > 




FIG. SI. Schematic diagram for control volume f2p and interpolation for a variable c based on the QUICK 
scheme. Here, e, to, n and s are the cell faces of the cell centered at P. 

are stored and small letter subscripts denotes the corresponding cell faces. Since the 
ion transport equations are coupled with the charge conservation equation, we adopt an 
iterative method. In order to reduce the system of algebraic equations into a tri-diagonal 
form, the variables at the locations 'EE' and 'W are taken as the previous iterated 
values. At every time-level, the solutions are obtained iteratively. 

At every time-level the iteration procedure starts with a guess for the electric potential 
at each cell center. At every iteration, the charge conservation equation is integrated over 
each control volume flp through the finite volume method. The elliptic PDE for charge 
conservation is solved by a line-by-line iterative method along with the successive-over- 
relaxation (SOR) technique. The iterations are continued until the absolute difference 
between two successive iterations becomes smaller than the tolerance limit 10 -6 for con- 
centration as well as for the electric potential. 
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S2 BOUNDARY CONDITIONS 




X 



FIG. S2. Schematic of the computational domain. 



The channel walls (boundaries 1 and 2 in figure S2 ) are impermeable for all ions 



N,.n = 0, V</>.n = 



(23) 



Here the subscript i — I, s, t refers to LE, sample and TE and n is the unit outward 
normal on the wall. 

Along the inlet (boundary 3) the concentrations and electric potential are prescribed 

as 



Hi + Ho Ht n 

Q = , Q = C s = U, 

Ht + Ho Hi 

while along the outlet (boundary 4) we set 

c t = c s = 0, Ct = 1, = 



E t H 



-x 3 



E,H 



-x A , 



(24) 



(25) 



where X3 = —X4 are the positions of the boundaries, located symmetrically around the 
origin. 



S3 EFFECT OF GRID SIZE, TIME EVOLUTION AND VALIDATION 



For code validation we consider pure ITP without a sample zone, i.e. only LE and TE 
present, since for this case an analytic result is available. This analytical solution due to 
Goet et al.^ for the TE concentration and axial electric field in a co-moving frame of 
reference under no bulk fluid flow can be expressed through the hypergeometric function 
as 

C t (x) 



C? 



111 + HQ AE X 

e ° 



and 



E t 



AE AE Ht + Ho 

C t (x) 



E(x) 



H1 + H0 m 
1 H e *o 
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(26) 



(27) 
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FIG. S3. Effect of grid size on the TE and LE concentration profiles and comparison with the analytical solution 
due to Goet et alP^. The electrophoretic mobility ratio is taken as ki(= /ii/fit = Di/D t ) = 2. 



At steady state the concentrations of TE and LE are related via 



where AE = E t — E\. Using relation (26), the LE concentration can easily be obtained 



from (28). 



Figure [S3] shows the effect of grid size on the TE and LE concentration profiles for a pure 
ITP without sample obtained with our code, together with the analytical results described 
above. As in the main text, the results are for Di = 7.0 x lCT 10 m 2 /s, E = 10 5 V/m, 
l t = li = 1/2 and a bulk LE concentration of Cf° = 10~ 3 M. The mobility ratio is taken 
to be k\ = Hi/fit — 2. Results are shown for three different grid sizes, with increasing 
fineness of grid in the X-direction. As can be seen, the analytical and numerical results 
agree very well for a grid of spacing SX = 5Y = 0.005. 



Figure S4 shows the time evolution towards the steady state for the case of ITP with 
a sample zone held stationary in the middle of the channel by a Poiseuille counterflow. 
As indicator the standard deviation, a, describing the width of the sample zone, is used 
here. The three cases shown correspond to values of the mobility ration k 2 = Hi/fi s close 
to the edges and in the middle of its range of validity, 1 < < k\. In all simulations 
shown in the main text care was taken that the steady state has indeed been reached. 

S4 TE, SAMPLE AND LE CONCENTRATION PROFILES 



For reference, we show in figure S5 the TE, sample and LE concentration profiles 
obtained for li = l t w 0.5 and l s ~ in the case without Poiseuille counterflow. The 
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FIG. S4. Time evolution of the standard deviation, a, reflecting the sample dispersion for ITP with Poiseuille 
counter-flow. Here fci(= Hi/fAt = Di/D t ) = 3 with k2(— ^i/fJ, s = Di/D s ) = 1.1, 2 and 2.9 and the amount of 
sample is C 3 /C ; °° = 40 fxm. 



cases shown range from sample amounts of C s /Cf° = 40 fim to C s /Cf° = Afim, two 
typical cases shown in the main text. In the former case the sample concentration clearly 
develops a plateau, while for the latter a Gaussian peak is observed. Note that both 
cases are in the regime of dispersed plateau mode, i.e. the area averaged concentration 
distribution shows a (distorted) Gaussian peak when Poiseuille counterflow is applied to 
keep the sample zone stationary. 

In the main text the primarily focus is on the distribution of the sample ions in the 
channel. An example for the effect on the LE and TE concentrations in the situation 



with Poiseuille counterflow is shown in figure S6 The area-averaged concentrations of all 
three ions (TE, sample and LE) is shown in the case of mobility ratios k\ = fii/nt — 3, 
^2 = fJ'i/fJ's — 2 and sample amount of C s /C^° = 40 /im for channel depths of H = 10 fim 
and 30 /im, respectively. As such the situation corresponds exactly to the one shown in 
figure 3 of the main text. It is evident, that the ID area-averaged model is able to capture 
the results obtained with the 2D model with a similar accuracy for the LE and TE as for 
the sample ions. 
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FIG. S5. Variation of TE (green), sample (red) and LE (blue) concentrations at steady state in a frame of 
reference moving at the constant speed U ITP for three different values of sample amount (a) Cs/C^ = 40 /im; 
(b) C s /Cf° — 10 fim and (c) C s /Cf° = 4/im. Here the mobility ratios are taken as fci(= fii/nt = Di/D t ) = 3 
and k 2 (= fli/fj, s = Di/D s ) = 2. 




(a) (b) 

FIG. S6. Area-averaged concentration distributions during ITP in Poiseuille counterflow of TE (green), sample 
(red) and LE (blue) for fci(= fii/lH = Di/Dt) = 3 and fe(= Mi/^s = Di/D a ) = 2 with channel depths of 
(a) H = 10 /im and (b) H = 30 /im. Full and dashed lines correspond to results obtained with the 2D and ID 
models, respectively. The amount of sample is C s /Cf° = 40 /im. (Cf. figure 3 in the main text, where only the 
sample distribution is shown.) 



